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1.  Introduction/Background 


CTH  is  a  continuum  hydrocode  originally  developed  by  Sandia  National  Laboratories.  The 
software  system  ean  model  multidimensional,  multimaterial,  large  deformation,  and  strong  shock 
wave  physics.^  The  existing  code  is  used  widely  across  the  US  Army  Research  Laboratory; 
increasing  the  code’s  capabilities  would  be  a  great  benefit.  For  example,  in  cases  that  involve 
material  fraeture,  continuum  Eulerian  codes  such  as  CTH  do  not  track  fractures  cleanly  and  can 
smear  eracks  and  change  their  extent  or  growth  eharacteristics.  These  adverse  numerieal  artifaets 
can  introduce  error  into  computational  results. 

The  reeent  implementation  of  the  material  point  method  (MPM)  into  CTH  may  ameliorate  some 
of  these  difficulties.  This  method  involves  the  discretization  of  solid  bodies  into  Lagrangian 
material  points.  These  points  (also  called  markers)  carry  the  full  state  of  the  material,  including 
mass,  density,  veloeity,  stress,  and  strain.  Interactions  among  adjacent  cells  of  the  numerical 
background  grid  enforce  compatibility  and  momentum  balance,  and  the  eonstitutive  relations  are 
evaluated  on  the  markers.  The  markers  receive  their  strain  rate  from  the  continuum  grid  during 
the  Lagrangian  part  of  the  step,  and  the  stress  and  material  history  are  updated  for  the  markers. 
The  material  points  remain  fixed  to  the  material  as  the  grid  is  returned  to  its  original  spatial 
loeation  during  a  remap  step.  The  remapped  markers  provide  the  stress  to  the  grid  for  the 
momentum  calculation  to  complete  the  eyele.  The  transient  eoupling  between  the  particles  and 
the  grid  enables  large  deformations  of  the  Lagrangian  particles  without  mesh  tangling.  Simply 
stated  here,  the  inclusion  of  MPM  has  the  potential  to  assist  in  aceurately  solving  all  kinds  of 
large  deformation  mechanics  problems  where  the  internal  state  of  material  is  history  dependent. 

The  incorporation  of  MPM  within  CTH  creates  a  new  eapability  to  be  evaluated  within  this 
report.  One  of  the  eentral  differences  between  continuum  CTH  and  CTH  with  MPM  is  the 
superposition  of  2  separate  pressure  terms  within  the  momentum  equation  (Eq.  1)  used  in  the 
analysis."^ 

*“  =  -V  •  (W  +  Q)dv  +  ^/ V  ■  (ff,  +  Pl)dV  +  ^fpsbdV.  (1) 

This  equation  can  be  broken  down  into  3  terms.  The  first  term  is  the  continuum  term;  it  is  the 
traditional  way  of  computing  accelerations  induced  by  the  pressure  gradient  on  the  background 
grid.  The  seeond  term  is  the  particle  term;  it  holds  the  full  stress  tensor,  at  the  particle  level. 
The  third  term  is  a  body  force  term,  whieh  does  not  play  a  role  in  the  following  computations. 
The  additive  contributions  of  the  positive  pressure  at  the  particle  level  and  the  negative  pressure 
at  the  continuum  level  play  a  significant  role  in  how  CTH  with  MPM  computes  aeeeleration  of 
the  grid.  In  a  uniform  field,  the  pressure  eontribution  removed  from  the  particle  term  is  exaetly 
replaced  by  the  continuum  pressure  from  the  baekground  grid.  Eor  inhomogeneous 
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deformation  fields,  the  result  is  more  eomplex,  as  will  be  suggested  in  this  work.  A  thorough 
understanding  of  the  CTH  solution  algorithm  with  markers  will  lead  to  greater  eonfidenee  for 
using  the  markers  in  a  wide  range  of  Army  problems. 


2.  Experiment/Calculations 


The  evaluation  of  CTH  with  MPM  is  eonducted  on  2  main  levels.  The  first  level  is  a  direct 
response  comparison  of  CTH  with  and  without  particles  for  a  zero-obliquity  impact  problem 
between  2  similar  materials.  The  second  level  of  the  evaluation  investigates  mesh  sensitivity  for 
the  particle  models  using  2  different  mesh  resolutions,  one  fine  and  another  coarse.  A 
comparison  between  these  2  resolutions  ensures  that  the  code  is  converging  to  a  physically 
reasonable  solution  based  on  different  levels  of  mesh  refinement.  Comparisons  between  the 
different  paired  simulations  briefly  described  here  demonstrate  how  accurately  the  new 
capabilities  of  CTH  operate.  The  input  deck  built  for  these  tests  is  composed  of  multiple 
components,  including  the  mesh,  diatom,  markers,  tracers,  and  the  equation  of  state  needed  to 
compute  the  stress  state  of  the  materials. 

2.1  Mesh 

The  mesh  is  a  user-supplied  spatial  discretization  of  the  material  geometries  for  the  problem  of 
interest.  For  a  Lagrangian  domain,  the  mesh  moves  with  the  material,  and  for  a  Eulerian  domain, 
the  material  advects  through  a  stationary  mesh.  Mesh  creation  for  the  current  evaluation  uses  a 
Eulerian  domain  with  a  2-dimensional  (2-D)  rectangular  geometry.  Its  mesh  extends  from  -0.01 
to  0.01  m  along  the  flyer  plate  propagation  direction.  There  are  400  elements  of  5  x  10'^  m 
dimension  for  the  fine  mesh  and  100  elements  of  2  x  10'"^  m  dimension  for  the  coarse  mesh.  The 
mesh  extends  0.002  m  from  the  origin  in  the  y-direction.  There  are  40  elements  of  5  x  10'^  m 
height  for  the  fine  mesh  and  10  elements  of  2  x  10'"^  m  height  for  the  coarse  mesh.  Symmetry 
conditions  were  applied  for  the  4  boundaries  of  the  computational  domain. 

2.2  Diatom 

The  diatom  block  of  the  input  deck  focuses  on  construction  of  the  bodies  of  material  that  make 
up  the  problem.  Here  the  initial  model  is  a  symmetric  2-D  plate  impact  experiment  between  2 
copper  plates.  The  plates  are  infinite  in  the  y-direction  and  measure  0.005  m  in  the  x-direction. 
They  travel  across  a  gap  of  0.0003  m  toward  each  other  at  500  m/s  each  until  impact. 

2.3  Particles 

The  continuum  solution  for  this  problem  is  evaluated  using  the  standard  CTH  code,  and  it  is  to 
be  compared  with  the  version  of  CTH  that  uses  the  MPM  (CTH  with  MPM)  particle  solution. 
The  particles  are  distributed  3  per  cell  in  the  x  and  y  directions  and  completely  represent  the 
copper  plates  described  in  the  diatom  for  the  particle  runs.  The  particles  are  capable  of 
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transporting  the  full  state  of  the  material  at  any  speeifie  partiele  loeation,  but  the  equation  of  state 
(EOS)  is  evaluated  on  the  background  grid.  The  particles  can  also  analyze  multiple  material 
pressures  and  temperatures  with  the  MMP4  option  enabled  in  this  portion  of  the  input  deck. 

2.4  Equation  of  State 

The  EOS  used  to  analyze  the  problem  is  the  Mie-Gruneisen  EOS.  The  strength  of  the  copper  is 
taken  from  the  Johnson-Cook  material  parameter  database  with  a  Poisson’s  ratio  of  0.3,  density 
of  8950  kg/m  ,  and  a  sound  speed  of  3560  m/s. 

2.5  Tracers 

Eight  tracers  are  originally  placed  along  the  lower  symmetry  boundary  of  the  left  copper  plate  in 
all  of  the  simulations.  These  tracers  point  to  the  locations  along  the  body  on  which  resultant  data 
will  be  saved,  recorded,  and  later  plotted.  Recorded  variables  on  the  tracers  include  time,  stress, 
position,  pressure,  temperature,  and  density.  Each  of  these  values  is  recorded  every  10  ns  during  the 
simulation. 

In  addition  to  using  an  appropriate  EOS  and  placing  tracers,  it  was  also  necessary  to  compute  the 
analytical  equilibrium  pressure  for  this  plate  impact.  This  calculation  establishes  a  standard 
baseline  around  which  computational  agreement  can  be  met.  This  pressure  was  evaluated  using 
Eq.  2.  It  multiplies  the  density,  p,  of  the  copper  plate  with  the  longitudinal  wave  speed.  Cl,  and 
velocity,  F,  of  the  plate. ^  An  equilibrium  pressure  of  2. 1  x  lO"^  MPa  was  calculated. 

P  =  P*Cl*V  .  (2) 

Not  only  should  the  numerically  calculated  peak  pressure  value  result  in  an  answer  on  the  order 
of  the  analytic  solution  (see  Eq.  2),  but  it  should  also  exhibit  a  waveform  like  the  one  shown  in 
Eig.  1 .  In  Eig.  I  a  continuum  impact  of  2  copper  plates  has  taken  place,  and  the  resulting  wave 
has  a  relatively  smooth  high-pressure  response  as  it  propagates  through  the  plates. 


Fig.  1  A  continuum  pressure  contour  plot  as  the  impact  wave  travels  after  plate  impact 
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3.  Results  and  Discussion 


3,1  Plate  Contact 

It  is  immediately  evident  that  CTH  with  MPM  is  vulnerable  to  one  of  the  inherent  problems  with 
eurrent  particle  methods;  void/particle  mixed  cell  issues.  Inspecting  the  pressure  contours  from 
CTH  for  the  initial  test  shows  that  the  2  copper  plates  are  not  entirely  coming  into  contact 
(Fig.  2).  However,  particles  away  from  the  impact  face  experience  a  pressure  increase  indicative 
of  impact.  The  mixed  cells  at  the  interface  of  the  2  plates  contain  both  particles  and  void  space. 
An  investigation  into  this  part  of  the  result  will  be  discussed  in  more  depth  for  a  less  complicated 
problem  later  in  this  report. 
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Fig.  2  Pressure  impact  image  from  CTH,  where  the  red  indicates  particle  areas  of  high  pressure, 
and  the  green  indicates  areas  of  zero  pressure  (scale  units  in  gigapascals) 


Upon  further  investigation,  it  was  also  discovered  that  the  particles  at  the  impact  face  were  at 
zero  pressure  despite  the  impact  having  already  taken  place,  as  shown  by  the  green  markers  in 
Fig.  2.  This  particular  issue  may  be  linked  to  the  superposition  of  the  particle  and  background 
grid  pressures  discussed  previously.  Determining  why  particles  are  indicating  zero  pressure  at 
the  impact  face  requires  further  investigation,  but  the  issue  associated  with  the  void  gap  between 
the  plates  is  mitigated  by  physically  closing  the  gap  within  the  problem  setup.  The  plates  are 
placed  directly  next  to  one  another  and  move  toward  each  other  for  impact. 

3,2  Lowered  Pressure  and  Mesh  Resolution 

With  the  void  gap  issues  having  been  noted  and  avoided,  it  is  now  possible  to  continue  the 
original  comparison  between  the  particle  method  of  CTH  with  MPM  and  continuum  CTH.  Even 
though  8  tracers  were  placed  along  the  left  plate  within  the  model,  only  tracer  4  will  be  used  for 
data  comparison  within  this  report  for  simplicity.  This  tracer  is  located  near  the  middle  of  the  left 
plate’s  horizontal  width,  and  it  is  indicated  by  the  yellow  circle  within  the  key  of  the  following 
figures.  The  trends  discussed  for  this  tracer  are  assumed  to  be  indicative  of  the  activity  at  the 
other  tracers. 
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The  left  graph  in  Fig.  3  depiets  the  pressure  impact  response  at  tracer  4  for  the  coarse  mesh 
compared  with  the  continuum  model.  For  this  calculation  the  response  of  the  particle  method  is 
well  below  that  of  the  continuum  solution  with  an  approximate  17%  difference  in  peak  pressure 
value.  It  is  also  important  to  note  that  the  continuum  solution’s  equilibrium  pressure  value  is 
approximately  equal  to  the  analytical  solution,  previously  calculated  as  2.08  x  10"^  MPa.  The  fine 
mesh  pressure  result  comparison  shown  in  the  right  graph  of  Fig.  3  has  better  agreement  with  the 
continuum  solution  with  only  a  half  percent  difference  in  peak  pressure;  however,  its  rise  time  is 
much  slower,  leading  to  a  pressure  difference  of  15%  after  the  initial  rise.  Direct  comparisons 
between  the  2  mesh  resolutions  for  CTH  with  MPM  are  shown  in  Fig.  4.  This  result  indicates 
that  the  solutions  from  these  2  simulations  are  not  identical.  There  is  approximately  a  14% 
difference  in  peak  pressure  values  between  the  2  resolution  models.  Additionally,  neither  of 
these  particle  results  exhibits  the  more  step-like  increase,  plateau,  and  step  decrease  in  pressure 
obtained  using  continuum  CTH. 
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Fig.  3  (Left)  Coarse  mesh  and  (right)  fine  mesh  continuum  comparisons 
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Fig.  4  Mesh  resolution  comparison  shows  disagreement 


3,3  Tracer  Location 

The  results  from  tracer  4  at  both  mesh  resolutions  are  generally  below  the  continuum  result.  In 
order  to  investigate  any  influence  from  the  boundary  conditions,  the  group  of  8  tracers  was 
shifted  vertically  up.  Now  instead  of  taking  resultant  data  from  the  lower  symmetry  boundary  of 
the  left  plate,  we  gather  it  from  the  middle  line  of  the  left  plate’s  geometry.  This  change  should 
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not  greatly  affect  the  resultant  data  from  the  tracers,  but  as  Fig.  5  shows,  the  shift  causes  a 
noticeable  increase  in  impact  pressure  by  approximately  7%.  This  indicates  that  there  are 
boundary  condition  effects  on  marker  tracers  placed  on  the  lower  symmetry  boundary.  This  is 
not  observed  using  continuum  CTH. 
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Fig.  5  Tracer  4  pressure  comparison  based  on  location 


3,4  Raised  Pressure 

The  change  in  pressure  output  for  tracer  4  as  a  result  of  its  change  in  position  requires  a 
reevaluation  of  the  original  comparison  between  CTH  with  MPM  and  continuum  CTH.  This  time 
we  focus  on  the  results  of  the  fine  mesh  since  it  more  closely  matches  the  continuum  results. 
Figure  6  shows  that  the  increase  in  pressure  response  obtained  by  tracer  relocation  has  placed  the 
particle  method  pressure  slightly  above  the  continuum  result  by  approximately  4%.  This  near 
reversal  of  the  original  pressure  discrepancy  can  be  qualified  further  only  by  investigating  the 
behavior  of  the  particle  method  specifically,  without  the  coupling  to  the  continuum  grid  through 
Eq.  1. 
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Fig.  6  Fine  mesh  continuum  comparison  after  the 
change  in  tracer  location 


3,5  Uintah:  A  Particle  Code 

To  gain  additional  insight  into  the  performance  of  the  CTH  with  MPM  particle  method,  a  second 
code  is  used  with  the  current  plate  impact  problem.  Uintah  is  a  pure  particle  code  developed  at 
the  University  of  Utah,  and  it  differs  from  CTH  in  that  it  does  not  factor  in  a  continuum  grid 
pressure  in  its  momentum  equation.^  The  original  symmetric  impact  across  a  gap  problem  is  built 


6 


within  Uintah’s  environment,  and  the  results  are  analyzed  to  elueidate  what  is  causing  the 
pressure  discrepancy  with  CTH  with  MPM. 

First,  the  issue  of  the  mixed  cells  halting  visible  plate  contact  was  examined.  This  behavior  can 
be  explained  by  the  method  by  which  particle  codes  move  bodies  like  the  copper  plates  within 
the  mesh.  Uintah  allows  the  plates  to  initially  move  toward  each  other  by  compressing  the  void 
space  in  between  them.  This  method  of  resolving  the  void  works  well  until  mixed  cells 
containing  both  particles  and  void  are  encountered  just  before  “physical”  impact.  Since  the 
particle  code  applies  the  same  deformation  rate  onto  copper  particles  and  void  within  any  cell, 
the  individual  particles  within  these  mixed  cells  are  compressed  just  like  void  space.  The  stress 
on  these  particles  increases  until  the  average  stress  on  their  cells  equals  the  equilibrium  pressure 
signaling  to  the  code  that  impact  has  now  occurred.  Essentially,  the  code  believes  that  the  plates 
have  come  in  contact  even  though  there  is  still  void  space  between  them. 

The  gap  issue  was  again  sidestepped  by  setting  the  2  plates  directly  next  to  each  other  for  the 
next  Uintah  simulation.  Fig.  7  shows  that  after  the  impact  of  the  2  plates,  the  resulting  waveform 
is  not  quite  smooth.  There  is  noticeable  ringing  occurring  as  the  wave  propagates  through  the 
plates.  A  more  in-depth  investigation  into  this  code’s  particle  behavior  is  necessary. 
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Fig.  7  Pressure  contour  from  Uintah  shows  ringing  as  the  impact  wave  travels  outward 

To  begin,  27  tracers  are  placed  at  consecutive  particle  locations  in  the  wake  of  the  impact  wave 
at  0.3  ps  after  impact.  The  pressure  response  of  these  tracers  shows  that  after  the  shock  passes,  2 
particles  in  a  given  cell  achieve  a  similar  stress  level  while  the  final  particle  falls  well  above  or 
below  in  stress.  It  seems  as  though  the  large  shock  caused  by  the  plate  impact  raises  the  stress  so 
quickly  with  respect  to  the  mesh  size  that  some  particles  comprising  the  plate  bodies  within  a  cell 
are  not  able  to  reach  the  final  equilibrium  pressure  state  before  being  pushed  into  an  adjoining 
cell.  However,  in  this  new  cell  the  existing  2  particles  have  already  achieved  their  final  state 
either  at  a  much  higher  or  lower  pressure.  The  result  is  a  set  of  high  and  low  pressures  that 
combines  to  create  the  correct  average  pressure.  As  this  effect  propagates  outward  from  the 
impact  face,  it  results  in  the  ringing  seen  in  the  pressure  plot  from  Uintah  in  Fig.  7. 

Placing  tracers  in  the  same  locations  for  the  same  time  step  in  CTH  with  MPM’s  input  deck 
yielded  similar  pressure  ringing  results  as  obtained  with  Uintah.  The  2  sets  of  data  were  then 
plotted  for  graphical  comparison.  The  scatter  plot  of  Fig.  8  shows  that  ringing  in  CTH  with  MPM 
is  still  apparent,  but  it  has  been  lessened  to  some  degree.  The  smoothing  of  the  CTH  with  MPM 
results  may  be  a  result  of  the  superposition  of  the  pressure  values  taken  from  the  individual 
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particles  and  the  baekground  grid  through  Eq.  1 .  This  mixing  of  particle  and  eontinuum  pressure 
is  not  done  in  Uintah.  Although  the  comparison  of  Uintah  as  a  particle  method  brought  about  an 
unexpeeted  result,  further  inspeetion  of  the  partiele  eode’s  output  brought  more  insight  into  what 
may  be  happening  within  the  partiele  implementation  of  CTH. 
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Fig.  8  Individual  particle  pressure  within  consecutive  cells  after  plate  impact  for  Uintah 
and  CTH  with  MPM 
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4.  Summary  and  Conclusions 


Several  issues  have  been  identified  within  the  new  implementation  of  the  material  point  method 
for  the  historically  continuum-based  hydrocode  CTH.  The  current  MPM  implementation  within 
CTH  exhibited  the  same  lack  of  mixed  cell  void  closure  between  impacted  bodies  as  Uintah. 
Boundary  condition  effects  on  particle  output  were  diseovered  upon  ehanging  tracer  locations 
within  the  model.  When  the  resolution  of  the  mesh  was  changed  from  coarse  to  fine,  the  pressure 
response  from  CTH  with  MPM  approached  the  CTH  continuum  solution.  This  mesh  sensitivity 
may  indieate  that  a  high  level  of  simulation  resolution  is  necessary  with  this  partiele 
implementation  in  order  to  achieve  the  stress  field  resolution  given  by  the  traditional  continuum 
method.  Also,  particles  appeared  to  exhibit  zero  pressure  values  at  the  impact  face  as  well  as 
noticeably  varied  values  within  a  single  eell  after  impact.  A  number  of  these  concerns  revolve 
around  the  addition  of  variables  evaluated  at  the  particle  and  background  grid  levels  within  CTH 
with  MPM’s  solution  algorithm.  Further  research  into  this  implementation  will  be  required  to 
gain  a  fundamental  understanding  of  how  this  averaging  method  eouples  to  partiele  models. 
Continued  development  of  CTH  with  MPM  will  hopefully  solve  the  problems  encountered 
during  the  application  of  this  2-D  plate  impact  problem  and  eventually  much  more  challenging 
meehanies  seenarios. 
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